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ABSTRACT 

This is the first of a series of papers aimed at developing and interpreting simulations 
of protoplanets interacting with turbulent accretion discs. In this hrst paper we study 
the turbulent disc models prior to the introduction of a perturbing protoplanet. We 
study cylindrical disc models in which a central domain is in Keplerian rotation and 
unstable to the magnetorotational instability (MRI). Models of varying disc size and 
aspect ratio H/r are considered with magnetic helds having zero net flux. We relate 
the properties of the turbulent models to classical viscous disc theory (Shakura & 
Sunyaev 1973). All models were found to attain a turbulent state in their Keplerian 
domains with volume averaged stress parameter a ^ 5 x 10“^. At any particular time 
the vertically and azimuthally averaged value exhibited large fluctuations in radius. 
However, an additional time average over periods exceeding 3 orbital periods at the 
outer boundary of the Keplerian domain resulted in a more smoothly varying quantity 
with radial variations within a factor of two or so. 

The vertically and azimuthally averaged radial velocity showed much larger spatial 
and temporal fluctuations, requiring additional time averaging for at least 7 — 8 orbital 
periods at the outer boundary of the Keplerian domain to limit them. Comparison with 
the value derived from the averaged stress using viscous disc theory yielded schematic 
agreement for feasible averaging times but with some indication that the effects of 
residual fluctuations remained. 

The behaviour described above must be borne in mind when considering laminar 
disc simulations with anomalous Navier-Stokes viscosity. This is because the operation 
of a viscosity as in classical viscous disc theory with anomalous viscosity coefficient 
cannot apply to a turbulent disc undergoing rapid changes due to external perturba¬ 
tion. The classical theory can only be used to describe the time averaged behaviour 
of the parts of the disc that are in a statistically steady condition for long enough for 
appropriate averaging to be carried out. 

Key words: accretion, accretion disks — MHD, instabilities, turbulence 


1 INTRODUCTION 

The recent and ongoing discovery of extrasolar giant plan¬ 
ets has stimulated renewed interest in the theory of planet 
formation (e.g. Mayor & Queloz 1995; Marcy, Cochran, & 
Mayor 1999; Vogt et al. 2002). In particular the discovery of 
giant planets close to their central stars has led to the idea 
that they migrated inwards due to gravitational interaction 
with the gaseous disc. Previous studies of the interaction 
between a protoplanet and a laminar but viscous disc [Pa¬ 
paloizou & Lin (1984); Lin & Papaloizou (1986, 1993); Bry- 
den et al. (1999); Nelson et al. (2000); D’Angelo, Henning & 
Kley (2002)] indicate that a protoplanet in the Jovian mass 
range will open a gap and that the torques exerted through 
the disc protoplanet interaction can produce inward orbital 


migration. However, the effect of the turbulence producing 
the anomalous viscosity has yet to be taken into account. 
The origin of this turbulence was uncertain until Balbus & 
Hawley (1991) provided an explanation for its origin through 
the operation of the magnetorotational instability (MRI). 
Improved computational resources now makes it feasible to 
consider three dimensional simulations of turbulent discs in¬ 
teracting with protoplanets. 

This is the first of a series of papers aimed at devel¬ 
oping and interpreting such simulations. In this first paper 
we focus on the turbulent disc models prior to the introduc¬ 
tion of a perturbing protoplanet. The effect of introducing 
the protoplanet will be considered in a companion paper 
(Nelson & Papaloizou (2002) - hereafter paper H). To ease 
computational requirements we adopt cylindrical disc mod- 
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els with no vertical stratification. For this first study we 
assume that the disc is adequately ionized throughout so 
that ideal MHD applies and consider models with zero net 
magnetic flux. We extend an initial study by Steinacker & 
Papaloizou (2002) (hereafter SP) to consider a wider range 
of models with varying disc size aspect ratio and resolution 
in the context of classical ‘q’ viscous disc theory (Shakura 
& Sunyaev 1973). This is because in addition to providing 
the most common and simple conceptual framework for disc 
modeling it still remains the main contact with observations 
(Balbus & Papaloizou 1999). We here focus on time averages 
of disc quantities such as the ref) stress and radial inflow ve¬ 
locity and study the extent to which stable behaviour of 
these quantities relates to classical disc theory. 

The plan of the paper is as follows. We describe the 
basic equations and mo del set up in §^. We describe thn 
numerical procedure in • 


!.l 


and our numerical results in 


Finally, we summarize our results in 


2 INITIAL MODEL SETUP 

The governing equations for MHD written in a frame rotat¬ 
ing with uniform angular velocity flpk with k being the unit 
vector in the vertical direction are: 

| + V.,v.O, (1) 

p + V • Vv^ + 2flpkxv = —Vp —pV'F 

+ ^(VxB)xB, (2) 

^=Vx(vxB). (3) 

where v, P, p, B and $ denote the fluid velocity, pressure, 

density, magnetic field, and potential, respectively. <E> con¬ 
tains contributions due to gravity and the centrifugal po¬ 
tential — {l/2)QpP. We use a locally isothermal equation of 
state 

P{r) = c{rf ■ p, (4) 

where c(r) denotes the sound speed which is specified as a 
fixed function of r. 

The models investigated may be described as cylindrical 
discs (e.g. Hawley 2001). Adopting cylindrical coordinates 
{z, r, <j)), the gravitational potential is taken to depend on r 
alone, so that 4> = —GM/r — (l/2)flpr^, where M is the 
central mass, G is the gravitational constant, and the sec¬ 
ond term represents the centrifugal potential where applica¬ 
ble. Thus the cylindrical disc models do not include a full 
treatment of the disc vertical structure. Models of this type 
are employed due to the high computational overhead that 
would be required to resolve fully the disc vertical structure 
of a stratified model. 

In this paper we consider the time dependent evolution 
and turbulent state that is set up in five models that we 
label as A, B, C, D, and E. These are all initiated with zero 
net magnetic flux in both vertical and azimuthal directions 
which is conserved for the duration of the simulations. As in 
SP periodic boundary conditions were used in the vertical 


and azimuthal directions, and each of the models has a ra¬ 
dial computational domain (ri,r 2 ) in which is embedded a 
central Keplerian domain (rai, ra 2 ) where the angular veloc¬ 
ity n oc and which becomes unstable due to the MRI. 

Interior and exterior to this Keplerian domain, adjacent to 
the inner and outer rigid radial boundaries, there exist re¬ 
gions in which the angular velocity profile is non Keplerian. 
These are described below. The relevant model parameters 
are given in table Dimensionless units of length and time 
are adopted such that ri = 1 and GM = 1. 

The angular velocity profile is chosen to be stable to 
the MRI in the boundary domains (where ri < r < rai and 
ra 2 < r < r 2 ), the inner one of which can be thought of 
as modeling the boundary layer between star and disc. For 
the models considered here, we adopted (?{r) cx 1/r with 
constant of proportionality such that c/(rD) =0.1 (models 
A, B, E) or c/{rQ,) = 0.2 (models C, D), apart from in 
the inner boundary domain where c^(r) oc In the 

Keplerian domain the initial density was such that p ocljr 
for all models. In models A, C, and D the angular velocity 
was initially constant in both boundary domains. In models 
B and E we took 12 oc in the inner boundary domain. 
Models A, C, and D were initiated with initial toroidal fields 
contained within the Keplerian domain while models B and 
E were initiated with poloidal fields. As in SP the initial 
magnetic field was given by 

Bi = Bo sin (2nR'K— —gj, (5) 

\ P ml Pm / 

with the integer Ur being the number of 27r cycles. The index 
i indicates either the vertical or the toroidal field component 
with the corresponding unit vector ej. Eor toroidal fields Bo 
is constant while for vertical fields Bo oc 1/r. The magnetic 
field was applied in an annulus within the Keplerian domain 
with inner and outer bounding radii Vm and r^i respectively. 
The normalization of Bq was chosen such that the initial 
magnetic energy in the Keplerian domain expressed in units 
of the volume integrated pressure there was 0.03 for models 
A, C, and D, 0.003 for model B, and 0.002 for model E. As 
the calculations proceed the magnetic field is seen to diffuse 
throughout the Keplerian domain until it is more or less 
fully magnetised, leading to a final turbulent state that is 
described in subsequent sections. 

Model E was initiated in the inertial frame with an az¬ 
imuthal domain of extent rr/S, and is the disc model used to 
study disc-planet interaction in paper II. This was run up 
to time t = 3825.5 when it had attained a fully turbulent 
state. Eor the purpose of studying the interaction between 
a turbulent protostellar disc and an embedded protoplanet, 
the azimuthal domain for this model was then extended to 
27r by stacking six of the n/3 sectors together. A transfor¬ 
mation to a rotating frame with Qp = 0.30645 was carried 
out and the evolution continued. Some results obtained at 
this stage are described below. We found no significant de¬ 
pendence of the evolution of the models on flp or the extent 
of the (p domain once the latter exceeded tt/3. 


2.1 Numerical procedure 

The numerical scheme that we employ is based on a spa¬ 
tially second-order accurate method that computes the ad- 
vection using the monotonic transport algorithm (Van Leer 
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Model 

Zl 

22 

ri 

r2 

02 

ra\ 

?*a2 

riz 

Ur 



‘^ml 

Ur 

A 

-0.2 

0.2 

1.0 

6.1 

7r/3 

1.25 

5.0 

40 

380 

100 

2.35 

4.35 

2 

B 

-0.2 

0.2 

1.0 

4.0 

n/2 

1.2 

3.7 

54 

334 

108 

1.33 

3.33 

6 

C 

-0.3 

0.3 

1.0 

4.5 

7r/3 

1.5 

3.7 

44 

334 

100 

2.33 

3.33 

1 

D 

-0.3 

0.3 

1.0 

4.5 

7r/3 

1.5 

3.7 

44 

334 

150 

2.33 

3.33 

1 

E 

-0.3 

0.3 

1.0 

00 

bo 

'K jZ 

1.2 

7.2 

60 

370 

100 

3.5 

6.5 

3 


Table 1. The first column gives the model label, the second and third give the vertical domain of extent the fourth and fifth give the 
radial domain while the sixth column gives the maximum extent of the azimuthal domain. The Keplerian domain is specified in columns 
seven and eight and the numbers of equally spaced grid points in the vertical, radial and azimuthal directions are given in columns nine, 
ten and eleven. The final three columns give the boundaries of the domain in which the initial magnetic field was applied and the number 
of 27r cycles in the functional form. 


1977). The MHD section of the code uses the method of 
characteristics constrained transport (MOCCT) as outlined 
in Hawley & Stone (1995) and implemented in the ZEUS 
code. The code has been developed from a version of NIR¬ 
VANA originally written by U. Ziegler (see SP, Ziegler & 
Rudiger (2000), and references therein) 


3 NUMERICAL RESULTS 

We now present numerical results for the evolution of models 
A - E into a saturated turbulent state and characterize their 
average evolutionary behaviour. To do this we use spatial 
and temporal averages of the state variables as indicated 
below. 


3.1 Vertically and horizontally averaged stresses 
and angular momentum transport 


In order to describe average properties of the turbulent 
models, we use quantities that are both vertically and az- 
imuthally averaged over the (()>, z) domain (e.g. Hawley 
2000) and in some cases an additional time average. The 
vertical and horizontal average of Q is defined through 


Q(r,t) 


f pQdzdcf) 
f pdzdfj) 


( 6 ) 


Below the average is taken over the full 27r in azimuth. 
If a smaller domain is used 27r should be replaced by the 
extent of this domain in what follows below in this section. 
In practice we have found that the results are independent 
of the size of the <!> domain, with the smallest azimuthal 
domain that we have considered being tt/S. 

The vertically and azimuthally averaged continuity equation 
(^) may then be written 

f + (7) 

at r or 

where the disc surface density is given by 


E = 


1 


pdzdcf). 


( 8 ) 


(Note that in performing the averaging the basic equations 
are integrated over the vertical and azimuthal domains with¬ 
out first multiplying by the density.) 


The vertically and azimuthally averaged Maxwell and 
Reynolds stresses, are respectively defined as follows: 

TM{r,t) = 27rEf Mz,r,<P,t)B,{z,rA,t) 

V 47rp 

and 

TRe{r, t) = 27vT,Svr{z, r, cj), t)Sv^{z, r, 0, t). (10) 

Here the velocity fluctuations Svr and Sv^ are defined 
through, 

Svr(z, r, <j), t) = Vr{z, T, <j), t) - tv(r, t), (11) 



5v^{z,r,(f),t) = v^{z,r,(j),t) - v^{r,t). 


( 12 ) 


The Shakura & Sunyaev (1973) a stress parameter appro¬ 
priate to the total stress is defined by 


a{r,t) 


Tfts — Tm 


(13) 


The vertically and azimuthally averaged azimuthal compo¬ 
nent of the equation of motion (2) can be written in the 
form 


d (Sj) 1 I d (r'Svrj) d (T,r^aP/p^ 
dt ^ r \ dr dr 


= 0 . 


(14) 


[see Balbus & Papaloizou (1999)]. Here j = rv^ is the spe¬ 
cific angular momentum. 

Using equation (|^, equation)]^) may also be written as 


Er 




d (Er^aP/p) 
dr 


(15) 


Note too that the averages may be extended without change 
of formalism to incorporate a time average such that for any 
quantity 

_ 1 _ 

<3(Di)^^y Q{r,t')dt', (16) 

where the time average is carried out over an interval 2A 
centered on the time t and this is incorporated into the def¬ 
inition of E. 

Under quasi-steady state conditions in the mean, if the 
time average is carried out over a sufficiently long interval, 
we should be able to neglect the time variation of the mean 
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Figure 1. Magnetic energy in the Keplerian domain expressed in 
units of the volume integrated pressure as a function of time for 
model A for t > 1378.6. This model was run for up to 2300 time 
units. This corresponds to 366 orbits at r = 1, 70 orbits at r = 3, 
and 33 orbits at r = 5. 


Specific angular momentum j. Then, as in the normal viscous 
disc theory, we expect to have an expression for the mean 
radial velocity of the form 


or or 


(17) 


3.2 Model A 

The time dependent evolution of the total magnetic energy 
in the Keplerian domain expressed in units of the volume 
integrated pressure, f^(B^/87r)dV/ PPdV, is plotted as a 
function of time for model A in figure [D. This model was run 
for up to 2300 time units. This corresponds to 366 orbits at 
r = 1, 70 orbits at r = 3, and 33 orbits at r = 5. The 
initial value of the ratio of the total magnetic energy to 
volume integrated pressure l/(/3) = J TP/{8'K)dV/ J PdV 
was 0.03 for model A and all others initiated with toroidal 
fields. After the onset of the MRI, some loss of the rather 
high initial magnetic energy due to reconnection occurs, and 
a relaxed turbulent state is attained after about ten orbits at 
the outer boundary of the Keplerian domain. The statistical 
properties of this do not depend on the initial conditions for 
models with zero net flux (see SP and below). 

The stress parameter a, defined in equation (^|), volume 
averaged over the Keplerian domain is plotted as a function 
of time for model A in figure Mean values are typically 
- 0.004. 

We now examine how the stress parameter a varies in 
space and time. To do this we consider time averages. Time 
averages of the stress parameter a are plotted as a function 
of dimensionless radius for model A in figure Time av¬ 
eraging ends at times t = 1391.6, t = 1430.0, t = 1549.7, 
t = 1646.4 and t = 1841.6, respectively. In each case the 
time averaging starts at t = 1378.6. Although a snapshot 
of a may reveal quite large variations (see figure below 
for model C, and also SP) after quite a short time a sta- 


1.01 


0 , 0 (« 



m 


(U-^^^^- 

m ISO! in m no * 

ME 

Figure 2. The stress parameter a volume averaged over the Ke¬ 
plerian domain is plotted as a function of time for model A. 


ble picture emerges. This is apparent after time averaging 
over an interval as short as 12 units, which represents 2 or¬ 
bits at the inner boundary and only 0.2 orbits at the outer 
boundary. Although from figure there is still some erratic 
behaviour visible for an averaging of 50 time units, for av¬ 
eraging periods exceeding 70 units, or one orbital period at 
the outer boundary of the active domain, the time averaged 
stress parameter appears to be a reasonably smooth func¬ 
tion of r. The variation of the time averaged a in the active 
domain is between 0.008 and 0.003. These values are typi¬ 
cal of those seen in local shearing box simulations starting 
with magnetic fields with zero net flux [Hawley, Gammie 
& Balbus (1996); Brandenburg et al. (1996); Brandenburg 
(1998); Fleming, Stone & Hawley (2000)]. Smaller values are 
obtained near the inner boundary. This is probably because 
magnetic field has not yet diffused to the inner boundary so 
that a long term turbulent steady state has not been reached 
there (note that the initial magnetic field was applied in an 
annulus away from the boundary domains as described in 
section 

In general the time averaged stress parameter a, and 
because of the relatively small fluctuations in time of E, 
also the time averaged stress itself attains a stable pattern 
after only a short period of time averaging. 

On the other hand, a stable pattern for the time aver¬ 
aged radial velocity takes significantly longer to attain. This 
is because from viscous disc theory we expect a characteris¬ 
tic value ~ 1.5a(H/r)^(rH) ~ 8 x A snapshot of 

the Vr averaged over the vertical and azimuthal domain is 
typically between one and two order of magnitudes higher, 
corresponding to large temporal fluctuations in that quan¬ 
tity (see also SP). 

We present a snapshot of = l/( 27 rLz) J pdfjjdz 

plotted as a function of radius for model A at the specific 
time t = 1841.6 in figure ^ (where Lz is the vertical extent 
of the computational domain). The temporal fluctuations in 
this quantity are found to be small, and representing the 
density averaged over the azimuthal and vertical domains it 
indicates some mass accretion towards the inner regions (see 
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Figure 3. Time averages of the stress parameter a are plotted 
as a function of dimensionless radius for model A. These end at 
times t = 1391.6 solid line, t = 1430.0 dashed line, t = 1549.7 
short dashed line t = 1646.4 dotted line and t = 1841.6 dot 
dashed line respectively. In each case the time averaging starts at 
t = 1378.6. 

also figure ^ for model E below). A snapshot of the product 
of the vertically and azimuthally averaged valnes of Vr with 
YjjLz, is plotted as a function of dimensionless radius for 
model A at time t = 1841.6 in figure This quantity is 
related to the instantaneous radial mass flux. Characteristic 
values of the averaged Vr are seen to be much larger than 
those expected from classical viscous disc theory which are 
recovered only after a long period of time averaging. 

To show this, we plot the product of the time averages of 
the vertically and azimuthally averaged values of Vr with the 
time averaged Tij L^, as a. function of dimensionless radius in 
figure ^ for model A. The averages end at times t = 1430.0, 
t = 1549.7, t = 1646.4 and t = 1841.6. The time averaging 
starts at t = 1378.6. Even the longest two averages over 463 
and 268 time units, although indicating magnitudes com¬ 
parable to those expected from viscous disc theory, deviate 
significantly. The shortest average over 52 time units is very 
different in character. From this we conclude that periods of 
time up to 7 orbital periods at the outer boundary of the ac¬ 
tive domain are required to obtain radial velocities that can 
be compared with viscous disc theory, and then the compar¬ 
ison can be moderately successful for the times over which 
averages can be taken here. 

We illustrate this by plotting the longest time average of 
the vertically and azimuthally averaged values of Vr with the 
time averaged T.jLz, as a function of dimensionless radius 
for model A in figure ^ The value of this quantity that we 
obtained from equation using the time averaged stress 



Figure 4. A snapshot of S/Lz is plotted as a function of radius 
for model A at the specific time t = 1841.6 



Figure 5. A snapshot of the product of the vertically and az¬ 
imuthally averaged values of Vr with HfLz, is plotted as a func¬ 
tion of dimensionless radius for model A at time t = 1841.6. 

Q.P/p obtained in the simulation, the Keplerian value for 
j, and using numerical differentiation is also plotted as the 
solid line. This latter quantity varies somewhat erratically 
because of the numerical differentiation, but nonetheless the 
general agreement is reasonable. 

3.3 Model B 

It is expected that models with initially zero net flux should 
attain similar steady states independent of initial conditions, 
disc thickness and size. To investigate this we considered 
model B which started from a vertical field and had outer 
boundary of the active domain at r = 3.7. The magnetic 
energy in the Keplerian domain expressed in units of the 
volume integrated pressure is plotted as a function of time 
for model B in figure |^. As for model A and all others to 
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Figure 6. The product of the time averages of the vertically 
and azimuthally averaged values of Vr with the time averaged 
Yi!Lz^ are plotted as a function of dimensionless radius for model 
A. These are taken at times t = 1430.0 (solid line), t = 1549.7 
(dashed line), t = 1646.4 (short dashed line) and t = 1841.6 
(dotted line). The time averaging starts from t = 1378.6. 



Figure 7. The product of the time averages of the vertically and 
azimuthally averaged values of Vr with the time averaged YjLz, 
is plotted as a function of dimensionless radius for model A in 
(dotted line). The time average starts at f = 1378.6 and ends 
at t = 1841.6. The value obtained from the time averaged stress 
using equation([l^ is also plotted (solid line). The method used 
for doing this is described in the text. 



0 100 200 300 400 500 600 700 800 

TIME 

Figure 8. Magnetic energy in the Keplerian domain expressed 
in units of the volume integrated pressure as a function of time 
for model B. 


0.01 
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Figure 9. The stress parameter a. volume averaged over the Ke¬ 
plerian domain is plotted as a function of time for model B. 


be presented, this settles down to a turbulent state with 
l/(/3) ~ 0.01. 

The stress parameter a, volume averaged over the Ke¬ 
plerian domain, is plotted as a function of time for model B 
in figure ^ This also takes on similar values and has similar 
behaviour to other models. 

The behaviour of the time averages of the stress pa¬ 
rameter a are plotted as a function of dimensionless radius 
for model B in figure |^. These end at times t = 525.4, 
t = 634.5, t = 664.3 and t = 740.6. The time averaging starts 
at t = 515.2. The shortest average over only 10 time units, 
gives erratic behaviour and even negative values. However, 
for those taken over more than 120 units a smooth stable 
behaviour is obtained with typical values of 0.004. But note 
that there are also long term cyclic trends which are seen in 
the time averages (see also Brandenburg 1998). 
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Figure 10. Time averages of the stress parameter a are plotted 
as a function of dimensionless radius for model B. These end at 
times t = 525.4 solid line, t = 634.5 dashed line, t = 664.3 short 
dashed line and t = 740.6 dotted line. The time averaging starts 
at t = 515.2. 


3.4 Models C and D 

We now describe models C and D. These models had an in¬ 
creased sound speed such that c/{rQ) — 0.2. Further, model 
D had 50 percent larger resolution in azimuth than model 
C in order to test the effects of azimuthal resolution. 

The magnetic energy in the Keplerian domain expressed 
in units of the volume integrated pressure is shown as a func¬ 
tion of time for models C and D in figure The behaviour 
of these models is very similar and again we find typical 
values ~ 0.01 for this quantity. 

Volume averages of the stress parameter a in the Kep¬ 
lerian domain for these models are plotted as a function of 
time in figure Again these models behave similarly. 

We plot time averages of the stress parameter a as a 
function of dimensionless radius for model C in figure [l3| 
These are taken at times t = 932.6, t = 983.0, t = 1024.8 
and t = 1163.4. The time averaging starts from t = 915.9. 
A similar plot for model D is given in figure |^. Here the 
averages are taken at times t = 833.6 solid line, and t = 
815.3 dashed line. The time averaging starts from t = 447.5. 

The time averages of a are more uniform in the thicker 
disc models C and D with stable behaviour being attained 
for averaging periods exceeding 4 orbital periods at the outer 
boundary of the Keplerian domain. The more uniform be¬ 
haviour may occur because the shorter viscous time in these 
cases enables a better relaxation of the disc to a state where 
memory of initial conditions is lost. 

To compare with the time averages a snapshot of the 
stress parameter a plotted as a function of radius for model 
C at time t = 1163.4 is shown in figure |^. 

To investigate the applicability of equation ( 0 ) , we plot 
the time average of the vertically and azimuthally averaged 


0.04 
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Figure 11. Magnetic energy in the Keplerian domain expressed 
in units of the volume integrated pressure as a function of time 
for model C solid line and model D dashed line. 



TIME 


Figure 12. Volume averages of the stress parameter a in the 
Keplerian domain as a function of time for model C solid line and 
model D dashed line. 

values of Vr with the time averaged as a function of 

dimensionless radius for model D in figure [l^ . The value 
of this quantity that we obtained from equation (0 , using 
the time averaged stress, the Keplerian value for j and us¬ 
ing numerical differentiation is also plotted as the solid line. 
As found for model A this quantity varies somewhat er¬ 
ratically because of the numerical differentiation, but there 
is schematic agreement. The time average was taken over 
eight orbital periods at the outer boundary of the Keplerian 
domain which is not enough to reduce the fluctuations in 
the radial velocity which can be up to two orders of mag¬ 
nitude larger than the mean to a very low level. This may 
require much longer averaging times comparable to the vis¬ 
cous timescale which are not feasible to consider here. 
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Figure 13. Time averages of the stress parameter a are plotted 
as a function of dimensionless radius for model C. These are taken 
at times t = 932.6 solid line, t = 983.0 dashed line, t = 1024.8 
short dashed line and t = 1163.4 dotted line. The time averaging 
starts from t = 915.9. 
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Figure 14. Time averages of the stress parameter a are plotted 
as a function of dimensionless radius for model D. These are taken 
at times t = 833.6 solid line, and t = 815.3 dashed line. The time 
averaging starts from t = 447.5. 


3.5 Model E 

This model was initiated with a poloidal magnetic field and 
had a value of c/{rQ) — 0.1 throughout the Keplerian do¬ 
main. The inner edge of the active Keplerian domain was 
located at r = 1.2 and the outer edge at r = 7.2, with the 
azimuthal extent being n/3, so this model has the largest 
radial extent that we consider in this work. This model is 
used in paper II to study the interaction between a MHD 
turbulent disc and a giant protoplanet. 

The temporal evolution of the magnetic energy in units 
of the volume integrated pressure is shown in figure This 
quantity saturates at a value of ^ 0.01 in the final turbulent 
state, similar to the previous runs described. The volume 
averaged stress parameter, a, is plotted as a function of time 
in figure ^ and shows similar behaviour to the previous 
runs described, saturating at a value of ~ 5 x 10“®. The 
radial variation of the time averaged stress parameter a is 
shown in figure with the solid fine denoting the total a, 
and the dashed fine denoting the magnetic contribution. The 
process of time averaging was initiated at a time t = 3825.5 
and completed at t = 4244.1. This time interval corresponds 
to 66.7 orbits at r = 1 and 3.5 orbits at r = 7.2. 

At a time of t = 3825.5, the disc model E was used 
to construct a model with an azimuthal extent of 27r for 
the purpose of studying disc-planet interactions, (whilst the 
original model E run was continued until t = 4500). This 
was achieved by simply patching six identical copies of the 
7r/3 model together to make a full 27r disc, ensuring of course 
that the azimuthal periodicity condition for the 7r/3 mod¬ 
els were properly utilised to construct the 27r disc. In ad¬ 
dition, the disc model was transformed from being in an 
inertial frame to a rotating frame with Qp = 0.30645. This 



Figure 15. A snapshot of the stress parameter a plotted as a 
function of radius for model C at time t = 1163.4 The solid curve 
corresponds to the total stress while the lower dashed curve is 
obtained when only the magnetic stresses are taken into account. 

value of Up corresponds to the angular velocity of material 
in circular Keplerian orbit at a radius of r = 2.2. We have 
checked that transforming the models from an inertial to a 
rotating frame, and/or increasing the azimuthal extent of a 
model by patching together identical copies, have no signif¬ 
icant effect on the statistical properties of the turbulence. 
Figures ^ and ^ show the magnetic energy (in units of 
the volume integrated pressure) and the volume integrated 
stress parameter a as functions of time for model E when 
transformed into the rotating frame (dotted line) and in the 
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Figure 16. The product of the time averages of the vertically 
and azimuthally averaged values of Vr with the time averaged 
TijLz, is plotted as a function of dimensionless radius for model 
D (dashed line). The time average starts at t = 447.5. and ends 
at t = 833.6 The value obtained from the time averaged stress 
using equation (|l7[) is also plotted (solid line). 


inertial frame (solid line). Small divergences in the numerical 
values are observed, since the code is now evolving a differ¬ 
ent numerical solution, but the statistical properties remain 
very similar. An almost identical situation arises when in¬ 
creasing the azimuthal domain from tt/S to 27r. 

A greyscale plot of the density at the vertical midplane 
of the disc is shown in figure corresponding to a time 
t = 3825.5. The usual trailing waves generated by the tur¬ 
bulence can clearly be seen in this figure. The variation of 
the azimuthally averaged midplane density as a function of 
radius is plotted in figure The solid line corresponds to 
the density distribution at time t = 3825.5, whereas the 
dashed line shows the initial values at t = 0. It is clear from 
this figure that some inward mass accretion has occurred. 


Magnetic Energy versus Time 



Figure 17. Magnetic energy in the Keplerian domain expressed 
in units of the volume integrated pressure as a function of time 
for model E. 


a versus Time 



Figure 18. The stress parameter a volume averaged over the 
Keplerian domain is plotted as a function of time for model E. 


4 DISCUSSION 

We have presented a study of cylindrical disc models in 
which a central domain in Keplerian rotation is unstable 
to the MRI. Models of varying disc size and aspect ra¬ 
tio H/r were considered. They were initiated with small 
scale magnetic fields with zero net flux. Conservation of 
poloidal and toroidal flux ensures that this situation is main¬ 
tained throughout the simulations. Input of flux through the 
boundaries which remain magnetically inactive does not oc¬ 
cur. The models all attain a turbulent state, with statistical 


properties that do not depend on the initial conditions which 
is expected for genuine dynamo action. 

As these models have been prepared for studies of disc- 
protoplanet interaction we have focused on relating the 
properties of the turbulent models to classical viscous disc 
theory [Shakura & Sunyaev (1973)]. This is an important 
issue, because besides providing a conceptual framework for 
disc studies, as emphasized by Balbus & Papaloizou (1999) 
it is still the main contact between disc theory and observa¬ 
tions. 

All models were found to attain a turbulent state with 
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Time Averoged a versus Rodlus 



Figure 19. The time averaged stress parameter a as a function 
of radius for the model E. 


Mognetic Energy versus Time 



3800 4000 4200 4400 4600 

Time 


Figure 20. Magnetic energy in the Keplerian domain expressed 
in units of the volume integrated pressure as a function of time 
for model E. The solid line denotes the calculation performed in 
the inertial frame, the dotted line denotes that calculated in the 
rotating frame. 

volume averaged stress parameter a ~ 5 x 10“® and mean 
/3“^ ~ 0.01. We also found that the same results were ob¬ 
tained in a rotating frame and in agreement with Hawley 
(2000) they were independent of the extent of the azimuthal 
domain when this exceeded tt/S. 

The vertically and azimuthally averaged stress param¬ 
eter showed large radial fluctuations. Time averaging for a 
period exceeding 3 orbital periods was found to signihcantly 
reduce them. For models with H/r = 0.1 stable variations 
with radius of a factor of two were then noted whereas for 
models with H/r = 0.2 less variation was seen. Variations 
in the time averaged quantities are most probably due to 
some memory of initial conditions which take up to a vis¬ 
cous time to eradicate, this being shorter for the thicker disc 
models. The higher resolution thicker disc model tended to 
have larger values of a than a lower resolution counterpart, 
which may reflect a residual dependence on resolution (see 
Brandenburg et al. 1996). 

The vertically and azimuthally averaged radial velocity 
showed large radial and temporal fluctuations of up to two 


a versus Time 



Figure 21. The volume averaged stress parameter a as a func¬ 
tion of time for the model E. The solid line denotes the calculation 
performed in the inertial frame, the dotted line denotes that cal¬ 
culated in the rotating frame. 


Surfoce Density of Turbulent Disc 



12 3 4 5 6 7 


Figure 22. This plot shows the variations in density at the ver¬ 
tical midplane the disc model E. 

orders of magnitude larger than the inflow velocity expected 
from classical viscous disc theory (see also SP). Time averag¬ 
ing for a period of at least 7 — 8 orbital periods at the outer 
boundary of the Keplerian domain was required to reveal 
values of a magnitude comparable to the expected viscous 
inflow velocity. Comparison with the value derived from the 
averaged stress using viscous disc theory yielded schematic 
agreement for feasible averaging times. It is likely that very 
long averaging times are needed to eliminate residual fluctu¬ 
ations in the mean radial velocity and that such an averaging 
operation may only be possible for a very quiet and thin disc 
that has relaxed to a long term statistical steady state. 

The behaviour described above must be borne in mind 
when considering laminar disc simulations with anoma- 
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Azimutholly Averoged Density versus Rodius 



Figure 23. This plot shows the radial variation of the az- 
imuthally averaged density at the vertical midplane of the disc 
model E. The solid line corresponds to time t = 3825.5, whereas 
the dashed line shows the initial values at i = 0. 

lous Navier-Stokes viscosity [e.g. Bryden et al (1999); Kley 
(1999); Lubow, Seibert & Artymowicz (1999); D’Angelo et al 
(2002)]. There, radial inflow velocities produced by the vis¬ 
cosity produce phenomena like mass flow through gaps when 
a protoplanet is embedded in a disc, due to the viscosity act¬ 
ing as a constantly acting source of friction. From the above 
discussion we might expect different dynamical behaviour in 
the gap region induced by a protoplanet in a genuinely tur¬ 
bulent disc (see paper II), where the instantaneous velocity 
fluctuations are much larger than the time averaged radial 
velocities arising from the turbulence-induced angular mo¬ 
mentum transport. More generally the classical viscous disc 
theory cannot apply to a disc undergoing rapid changes due 
to external perturbation. It can only be used to describe the 
parts of the disc that are in a quasi-steady condition for long 
enough for appropriate averaging to be carried out. 
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